!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>

!     the subroutine LIMBES_feq computes the equilibrium distribution function feq
!     depends of the order Hermite expansion LIMBES_NH
!     returned argument
!     feq(1:LIMBES_nv) -> equilibrium distribution function (real*8)
!     needed arguments
!     real*8 rho macroscopic density
!     real*8 ux,uy macroscopic velocity component
!     real*8 T macroscopic temperature
subroutine LIMBES_calc_feq(rho,uux,uuy,RT,feq)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none

  integer*8 :: vi,NHx,NHy
  real*8 :: rho,ux,uy,uux,uuy,RT,teta
  real*8 :: vdotv,vdotu,udotu
  real*8 :: feq0,vx,vy,aa,HH,fact
  real*8 :: feq(1:LIMBES_nv)

  teta=RT/FLUID_RT0
  udotu=(uux*uux + uuy*uuy)/FLUID_RT0
  ux=uux/dsqrt(FLUID_RT0)
  uy=uuy/dsqrt(FLUID_RT0)

  do vi=1,LIMBES_nv
     feq0=LIMBES_feq0(vi)
     vx=LIMBES_vx(vi)/dsqrt(FLUID_RT0)
     vy=LIMBES_vy(vi)/dsqrt(FLUID_RT0)
     vdotv=(vx*vx + vy*vy)/FLUID_RT0
     vdotu=(vx*ux + vy*uy)/FLUID_RT0
     feq(vi)=zero
     do NHx=0,LIMBES_NH
        do NHy=0,LIMBES_NH
           feq(vi)= feq(vi) + (one/(fact(NHx)*fact(NHy)))*aa(NHx,ux,teta)*aa(NHy,uy,teta)*HH(NHx,vx)*HH(NHy,vy)   
        enddo
     enddo
     feq(vi)=rho*feq0*feq(vi)
  enddo

end subroutine LIMBES_calc_feq

!     ========================================================
!     ========================================================
!     ========================================================     
real*8 function feq_vi(vi,rho,uux,uuy,RT)
  use LIMBES_mod_fluid_config
  use LIMBES_mod_var
  implicit none

  integer*8 :: vi,NHx,NHy
  real*8 :: rho,ux,uy,uux,uuy,RT,teta
  real*8 :: vdotv,vdotu,udotu
  real*8 :: feq0,vx,vy,aa,HH,fact,temp_feq

  teta=RT/FLUID_RT0
  udotu=(uux*uux + uuy*uuy)/FLUID_RT0
  ux=uux/dsqrt(FLUID_RT0)
  uy=uuy/dsqrt(FLUID_RT0)


  feq0=LIMBES_feq0(vi)
  vx=LIMBES_vx(vi)/dsqrt(FLUID_RT0)
  vy=LIMBES_vy(vi)/dsqrt(FLUID_RT0)
  vdotv=(vx*vx + vy*vy)/FLUID_RT0
  vdotu=(vx*ux + vy*uy)/FLUID_RT0
  temp_feq=zero
  do NHx=0,LIMBES_NH
     do NHy=0,LIMBES_NH
        temp_feq= temp_feq + (one/(fact(NHx)*fact(NHy)))*aa(NHx,ux,teta)*aa(NHy,uy,teta)*HH(NHx,vx)*HH(NHy,vy)   
     enddo
  enddo
  temp_feq=rho*feq0*temp_feq
  feq_vi=temp_feq

end function feq_vi

!     ========================================================
!     ========================================================
!     ========================================================  
real*8 FUNCTION FACT(N)
  implicit none
  integer*8 ::  N,J
  FACT=1.0d+0
  DO J=2,N
     FACT=FACT*dble(J)
  enddo
END FUNCTION FACT

!     ========================================================
!     ========================================================
!     ========================================================  
!!$
!!$real*8 function HH(q,v)
!!$  implicit none
!!$  integer*8 :: q
!!$  real*8 :: v
!!$
!!$  select case(q)
!!$  case(0)
!!$     HH=1.0d+0
!!$
!!$  case(1)
!!$     HH=v
!!$
!!$  case(2)
!!$     HH=v**2-1
!!$
!!$  case(3)
!!$     HH=v**3-3*v
!!$
!!$  case(4)
!!$     HH=v**4-6*v**2+3
!!$
!!$  case(5)
!!$     HH=v**5-10*v**3+15*v
!!$
!!$  case(6)
!!$     HH=v**6-15*v**4+45*v**2-15
!!$
!!$  case(7)
!!$     HH=v**7-21*v**5+105*v**3-105*v
!!$
!!$  case(8)
!!$     HH=v**8-28*v**6+210*v**4-420*v**2+105  
!!$
!!$  case(9)
!!$     HH=v**9-36*v**7+378*v**5-1260*v**3+945*v
!!$
!!$  case(10)
!!$     HH=v**10-45*v**8+630*v**6-3150*v**4+4725*v**2-945
!!$
!!$  case(11)
!!$     HH=v**11-55*v**9+990*v**7-6930*v**5+17325*v**3-10395*v
!!$
!!$  case(12)
!!$     HH=v**12-66*v**10+1485*v**8-13860*v**6+51975*v**4-62370*v**2 +10395
!!$
!!$  case(13)
!!$     HH= v**13-78*v**11+2145*v**9-25740*v**7+135135*v**5-270270*v**3+135135*v
!!$
!!$  case(14)
!!$     HH= v**14-91*v**12+3003*v**10-45045*v**8+315315*v**6-945945*v**4+945945*v**2-135135 
!!$
!!$  case(15)
!!$     HH= v**15-105*v**13+4095*v**11-75075*v**9+675675*v**7-2837835*v**5+4729725*v**3-2027025*v
!!$
!!$  case(16)
!!$     HH= v**16-120*v**14+5460*v**12-120120*v**10+1351350*v**8-7567560*v**6+18918900*v**4-16216200*v**2+2027025
!!$
!!$  case(17)
!!$     HH= v**17-136*v**15+7140*v**13-185640*v**11+2552550*v**9-18378360*v**7+64324260*v**5-91891800*v**3+34459425*v
!!$
!!$  case(18)
!!$     HH= v**18-153*v**16+9180*v**14-278460*v**12+4594590*v**10-41351310*v**8+192972780*v**6-413513100*v**4+310134825*v**2-34459425
!!$
!!$  case(19)
!!$     HH= v**19-171*v**17+11628*v**15-406980*v**13+7936110*v**11-87297210*v**9+523783260*v**7-1571349780*v**5+1964187225*v**3-654729075*v
!!$
!!$  case(20)
!!$     HH= v**20-190*v**18+14535*v**16-581400*v**14+13226850*v**12-174594420*v**10+1309458150*v**8-5237832600*v**6+9820936125*v**4-6547290750*v**2+654729075
!!$  end select
!!$
!!$end function HH
!!$
!!$!     -------------------------------------------------------------------
!!$!     -------------------------------------------------------------------
!!$!     -------------------------------------------------------------------
!!$
!!$real*8 function aa(q,u,T)
!!$  implicit none
!!$  integer*8 :: q
!!$  real*8 :: u,T
!!$
!!$  select case(q)
!!$  case(0)
!!$     aa=1.0d+0
!!$
!!$  case(1)
!!$     aa=u
!!$
!!$  case(2)
!!$     aa=T+u**2-1
!!$
!!$  case(3)
!!$     aa=3*u*T+u**3-3*u
!!$
!!$  case(4)
!!$     aa=3*T**2+(6*u**2-6)*T+u**4-6*u**2+3
!!$
!!$  case(5)
!!$     aa=15*u*T**2+(10*u**3-30*u)*T+u**5-10*u**3+15*u
!!$
!!$  case(6)
!!$     aa= 15*T**3+(45*u**2-45)*T**2+(15*u**4-90*u**2+45)*T+u**6-15*u**4+45*u**2-15
!!$
!!$  case(7)
!!$     aa= 105*u*T**3+(105*u**3-315*u)*T**2+(21*u**5-210*u**3+315*u)*T+u**7-21*u**5+105*u**3-105*u
!!$
!!$  case(8)
!!$     aa= 105*T**4+(420*u**2-420)*T**3+(210*u**4-1260*u**2+630)*T**2+(28*u**6-420*u**4+1260*u**2-420)*T+u**8-28*u**6+210*u**4-420*u**2+105
!!$
!!$  case(9)
!!$     aa= 945*u*T**4+(1260*u**3-3780*u)*T**3+(378*u**5-3780*u**3+5670*u)*T**2+(36*u**7-756*u**5+3780*u**3-3780*u)*T+u**9-36*u**7+378*u**5-1260*u**3+945*u      
!!$
!!$  case(10)
!!$     aa= 945*T**5+(4725*u**2-4725)*T**4+(3150*u**4-18900*u**2+9450)*T**3+(630*u**6-9450*u**4+28350*u**2-9450)*T**2+(45*u**8-1260*u**6+9450*u**4-18900*u**2+4725)*T+u**10-45*u**8+630*u**6-3150*u**4+4725*u**2-945
!!$
!!$  case(11)
!!$     aa= 10395*u*T**5+(17325*u**3-51975*u)*T**4+(6930*u**5-69300*u**3+103950*u)*T**3+(990*u**7-20790*u**5+103950*u**3-103950*u)*T**2+(55*u**9-1980*u**7+20790*u**5-69300*u**3+51975*u)*T+u**11-55*u**9+990*u**7-6930*u**5+17325*u**3-10395*u      
!!$
!!$  case(12)
!!$     aa= 10395*T**6+(62370*u**2-62370)*T**5+(51975*u**4-311850*u**2+155925)*T**4+(13860*u**6-207900*u**4+623700*u**2-207900)*T**3+(1485*u**8-41580*u**6+311850*u**4-623700*u**2+155925)*T**2+(66*u**10-2970*u**8+41580*u**6-207900*u**4+311850*u**2-62370)*T+u**12-66*u**10+1485*u**8-13860*u**6+51975*u**4-62370*u**2+10395
!!$
!!$  case(13)
!!$     aa= 135135*u*T**6+(270270*u**3-810810*u)*T**5+(135135*u**5-1351350*u**3+2027025*u)*T**4+(25740*u**7-540540*u**5+2702700*u**3-2702700*u)*T**3+(2145*u**9-77220*u**7+810810*u**5-2702700*u**3+2027025*u)*T**2+(78*u**11-4290*u**9+77220*u**7-540540*u**5+1351350*u**3-810810*u)*T+u**13-78*u**11+2145*u**9-25740*u**7+135135*u**5-270270*u**3+135135*u
!!$
!!$  case(14)
!!$     aa= 135135*T**7+(945945*u**2-945945)*T**6+(945945*u**4-5675670*u**2+2837835)*T**5+(315315*u**6-4729725*u**4+14189175*u**2-4729725)*T**4+(45045*u**8-1261260*u**6+9459450*u**4-18918900*u**2 +4729725)*T**3+(3003*u**10-135135*u**8+1891890*u**6-9459450*u**4+14189175*u**2-2837835)*T**2+(91*u**12-6006*u**10+135135*u**8-1261260*u**6+4729725*u**4-5675670*u**2+945945)*T+u**14-91*u**12+3003*u**10-45045*u**8+315315*u**6-945945*u**4+945945*u**2-135135
!!$
!!$  case(15)
!!$     aa= 2027025*u*T**7+(4729725*u**3-14189175*u)*T**6+(2837835*u**5-28378350*u**3+42567525*u)*T**5+(675675*u**7-14189175*u**5+70945875*u**3-70945875*u)*T**4+(75075*u**9-2702700*u**7+28378350*u**5-94594500*u**3+70945875*u)*T**3+(4095*u**11-225225*u**9+4054050*u**7-28378350*u**5+70945875*u**3-42567525*u)*T**2+(105*u**13-8190*u**11+225225*u**9-2702700*u**7+14189175*u**5-28378350*u**3+14189175*u)*T+u**15-105*u**13+4095*u**11-75075*u**9+675675*u**7-2837835*u**5+4729725*u**3-2027025*u 
!!$
!!$  case(16)
!!$     aa= 2027025*T**8+(16216200*u**2-16216200)*T**7+(18918900*u**4-113513400*u**2+56756700)*T**6+(7567560*u**6-113513400*u**4+340540200*u**2-113513400)*T**5+(1351350*u**8-37837800*u**6+283783500*u**4-567567000*u**2+141891750)*T**4+(120120*u**10-5405400*u**8+75675600*u**6-378378000*u**4+567567000*u**2-113513400)*T**3+(5460*u**12-360360*u**10+8108100*u**8-75675600*u**6+283783500*u**4-340540200*u**2+56756700)*T**2+(120*u**14-10920*u**12+360360*u**10-5405400*u**8+37837800*u**6-113513400*u**4+113513400*u**2 -16216200)*T+u**16-120*u**14+5460*u**12-120120*u**10+1351350*u**8-7567560*u**6+18918900*u**4-16216200*u**2+2027025
!!$
!!$  case(17)
!!$     aa= 34459425*u*T**8+(91891800*u**3-275675400*u)*T**7+(64324260*u**5-643242600*u**3+964863900*u)*T**6+(18378360*u**7-385945560*u**5+1929727800*u**3-1929727800*u)*T**5+(2552550*u**9-91891800*u**7+964863900*u**5-3216213000*u**3+2412159750*u)*T**4+(185640*u**11-10210200*u**9+183783600*u**7-1286485200*u**5+3216213000*u**3-1929727800*u)*T**3+(7140*u**13-556920*u**11+15315300*u**9-183783600*u**7+964863900*u**5-1929727800*u**3+964863900*u)*T**2+(136*u**15-14280*u**13+556920*u**11-10210200*u**9+91891800*u**7-385945560*u**5+643242600*u**3-275675400*u)*T+u**17-136*u**15 +7140*u**13-185640*u**11+2552550*u**9-18378360*u**7+64324260*u**5-91891800*u**3+34459425*u
!!$
!!$  case(18)
!!$     aa= 34459425*T**9+(310134825*u**2-310134825)*T**8+(413513100*u**4-2481078600*u**2+1240539300)*T**7+(192972780*u**6-2894591700*u**4+8683775100*u**2-2894591700)*T**6+(41351310*u**8-1157836680*u**6+8683775100*u**4-17367550200*u**2+4341887550)*T**5+(4594590*u**10-206756550*u**8+2894591700*u**6-14472958500*u**4+21709437750*u**2-4341887550)*T**4+(278460*u**12-18378360*u**10+413513100*u**8-3859455600*u**6+14472958500*u**4-17367550200*u**2+2894591700)*T**3+(9180*u**14-835380*u**12+27567540*u**10-413513100*u**8+2894591700*u**6-8683775100*u**4+8683775100*u**2-1240539300)*T**2+(153*u**16-18360*u**14+835380*u**12-18378360*u**10+206756550*u**8-1157836680*u**6+2894591700*u**4-2481078600*u**2+310134825)*T+u**18-153*u**16+9180*u**14-278460*u**12+4594590*u**10- 41351310*u**8+192972780*u**6-413513100*u**4+310134825*u**2-34459425
!!$
!!$  case(19)
!!$     aa= 654729075*u*T**9+(1964187225*u**3-5892561675*u)*T**8+(1571349780*u**5-15713497800*u**3+23570246700*u)*T**7+(523783260*u**7-10999448460*u**5+54997242300*u**3-54997242300*u)*T**6+(87297210*u**9-3142699560*u**7+32998345380*u**5-109994484600*u**3+82495863450*u)*T**5+(7936110*u**11-436486050*u**9+7856748900*u**7-  54997242300*u**5+137493105750*u**3-82495863450*u)*T**4+(406980*u**13-31744440*u**11+872972100*u**9-10475665200*u**7+54997242300*u**5-109994484600*u**3+54997242300*u)*T**3+(11628*u**15-1220940*u**13+47616660*u**11-872972100*u**9+7856748900*u**7-32998345380*u**5+54997242300*u**3-23570246700*u)*T**2+(171*u**17-23256*u**15+1220940*u**13-31744440*u**11+436486050*u**9-3142699560*u**7+10999448460*u**5-15713497800*u**3+5892561675*u)*T+u**19-171*u**17+11628*u**15-406980*u**13+7936110*u**11-87297210*u**9+523783260*u**7-1571349780*u**5+1964187225*u**3-654729075*u
!!$
!!$  case(20)
!!$     aa= 654729075*T**10+(6547290750*u**2-6547290750)*T**9+(9820936125*u**4-58925616750*u**2+29462808375)*T**8+(5237832600*u**6-78567489000*u**4+235702467000*u**2-78567489000)*T**7+(1309458150*u**8-36664828200*u**6+274986211500*u**4-549972423000*u**2+137493105750)*T**6+(174594420*u**10-7856748900*u**8+109994484600*u**6-549972423000*u**4+824958634500*u**2-164991726900)*T**5+(13226850*u**12-872972100*u**10+19641872250*u**8-183324141000*u**6+687465528750*u**4-824958634500*u**2+137493105750)*T**4+(581400*u**14-52907400*u**12+1745944200*u**10-26189163000*u**8+183324141000*u**6-549972423000*u**4+549972423000*u**2-78567489000)*T**3+ (14535*u**16-1744200*u**14+79361100*u**12-1745944200*u**10+19641872250*u**8-109994484600*u**6+274986211500*u**4-235702467000*u**2+29462808375)*T**2+(190*u**18-29070*u**16+1744200*u**14-52907400*u**12+872972100*u**10-7856748900*u**8+36664828200*u**6-78567489000*u**4+58925616750*u**2-6547290750)*T+u**20-190*u**18+14535*u**16-581400*u**14+13226850*u**12-174594420*u**10+1309458150*u**8-5237832600*u**6+9820936125*u**4-6547290750*u**2+654729075     
!!$
!!$  end select
!!$end function aa

!     ========================================================
!     ========================================================
!     ========================================================  

real*8 function HH(q,v)
  implicit none
  integer*8 :: q
  real*8 :: v


  if (q.EQ.0) then
     HH=1.0d+0
  elseif(q.EQ.1) then
     HH=v
  elseif(q.EQ.2) then
     HH=v**2-1
  elseif(q.EQ.3) then 
     HH=v**3-3*v
  elseif(q.EQ.4) then
     HH=v**4-6*v**2+3
  elseif(q.EQ.5) then
     HH=v**5-10*v**3+15*v
  elseif(q.EQ.6) then
     HH=v**6-15*v**4+45*v**2-15
  elseif(q.EQ.7) then
     HH=v**7-21*v**5+105*v**3-105*v
  elseif(q.EQ.8) then
     HH=v**8-28*v**6+210*v**4-420*v**2+105  
  elseif(q.EQ.9) then
     HH=v**9-36*v**7+378*v**5-1260*v**3+945*v
  elseif(q.EQ.10) then
     HH=v**10-45*v**8+630*v**6-3150*v**4+4725*v**2-945
  elseif(q.EQ.11) then
     HH=v**11-55*v**9+990*v**7-6930*v**5+17325*v**3-10395*v
  elseif(q.EQ.12) then
     HH=v**12-66*v**10+1485*v**8-13860*v**6+51975*v**4-62370*v**2 +10395
  elseif(q.EQ.13) then
     HH= v**13-78*v**11+2145*v**9-25740*v**7+135135*v**5-270270*v**3+135135*v
  elseif(q.EQ.14) then
     HH= v**14-91*v**12+3003*v**10-45045*v**8+315315*v**6-945945*v**4+945945*v**2-135135 
  elseif(q.EQ.15) then
     HH= v**15-105*v**13+4095*v**11-75075*v**9+675675*v**7-2837835*v**5+4729725*v**3-2027025*v
  elseif(q.EQ.16) then
     HH= v**16-120*v**14+5460*v**12-120120*v**10+1351350*v**8-7567560*v**6+18918900*v**4-16216200*v**2+2027025
  elseif(q.EQ.17) then
     HH= v**17-136*v**15+7140*v**13-185640*v**11+2552550*v**9-18378360*v**7+64324260*v**5-91891800*v**3+34459425*v
  elseif(q.EQ.18) then
     HH= v**18-153*v**16+9180*v**14-278460*v**12+4594590*v**10-41351310*v**8+192972780*v**6-413513100*v**4+310134825*v**2-34459425
  elseif(q.EQ.19) then
     HH= v**19-171*v**17+11628*v**15-406980*v**13+7936110*v**11-87297210*v**9+523783260*v**7-1571349780*v**5+1964187225*v**3-654729075*v
  elseif(q.EQ.20) then
     HH= v**20-190*v**18+14535*v**16-581400*v**14+13226850*v**12-174594420*v**10+1309458150*v**8-5237832600*v**6+9820936125*v**4-6547290750*v**2+654729075
  endif

end function HH

!     -------------------------------------------------------------------
!     -------------------------------------------------------------------
!     -------------------------------------------------------------------

real*8 function aa(q,u,T)
  implicit none
  integer*8 :: q
  real*8 :: u,T

  if (q.EQ.0) then
     aa=1.0d+0
  elseif (q.EQ.1) then
     aa=u
  elseif (q.EQ.2) then
     aa=T+u**2-1
  elseif (q.EQ.3) then
     aa=3*u*T+u**3-3*u
  elseif (q.EQ.4) then
     aa=3*T**2+(6*u**2-6)*T+u**4-6*u**2+3
  elseif (q.EQ.5) then
     aa=15*u*T**2+(10*u**3-30*u)*T+u**5-10*u**3+15*u
  elseif (q.EQ.6) then
     aa= 15*T**3+(45*u**2-45)*T**2+(15*u**4-90*u**2+45)*T+u**6-15*u**4+45*u**2-15
  elseif (q.EQ.7) then
     aa= 105*u*T**3+(105*u**3-315*u)*T**2+(21*u**5-210*u**3+315*u)*T+u**7-21*u**5+105*u**3-105*u
  elseif (q.EQ.8) then
     aa= 105*T**4+(420*u**2-420)*T**3+(210*u**4-1260*u**2+630)*T**2+(28*u**6-420*u**4+1260*u**2-420)*T+u**8-28*u**6+210*u**4-420*u**2+105
  elseif (q.EQ.9) then
     aa= 945*u*T**4+(1260*u**3-3780*u)*T**3+(378*u**5-3780*u**3+5670*u)*T**2+(36*u**7-756*u**5+3780*u**3-3780*u)*T+u**9-36*u**7+378*u**5-1260*u**3+945*u     
  elseif (q.EQ.10) then
     aa= 945*T**5+(4725*u**2-4725)*T**4+(3150*u**4-18900*u**2+9450)*T**3+(630*u**6-9450*u**4+28350*u**2-9450)*T**2+(45*u**8-1260*u**6+9450*u**4-18900*u**2+4725)*T+u**10-45*u**8+630*u**6-3150*u**4+4725*u**2-945
  elseif (q.EQ.11) then
     aa= 10395*u*T**5+(17325*u**3-51975*u)*T**4+(6930*u**5-69300*u**3+103950*u)*T**3+(990*u**7-20790*u**5+103950*u**3-103950*u)*T**2+(55*u**9-1980*u**7+20790*u**5-69300*u**3+51975*u)*T+u**11-55*u**9+990*u**7-6930*u**5+17325*u**3-10395*u    
  elseif (q.EQ.12) then
     aa= 10395*T**6+(62370*u**2-62370)*T**5+(51975*u**4-311850*u**2+155925)*T**4+(13860*u**6-207900*u**4+623700*u**2-207900)*T**3+(1485*u**8-41580*u**6+311850*u**4-623700*u**2+155925)*T**2+(66*u**10-2970*u**8+41580*u**6-207900*u**4+311850*u**2-62370)*T+u**12-66*u**10+1485*u**8-13860*u**6+51975*u**4-62370*u**2+10395
  elseif (q.EQ.13) then
     aa= 135135*u*T**6+(270270*u**3-810810*u)*T**5+(135135*u**5-1351350*u**3+2027025*u)*T**4+(25740*u**7-540540*u**5+2702700*u**3-2702700*u)*T**3+(2145*u**9-77220*u**7+810810*u**5-2702700*u**3+2027025*u)*T**2+(78*u**11-4290*u**9+77220*u**7-540540*u**5+1351350*u**3-810810*u)*T+u**13-78*u**11+2145*u**9-25740*u**7+135135*u**5-270270*u**3+135135*u
  elseif (q.EQ.14) then
     aa= 135135*T**7+(945945*u**2-945945)*T**6+(945945*u**4-5675670*u**2+2837835)*T**5+(315315*u**6-4729725*u**4+14189175*u**2-4729725)*T**4+(45045*u**8-1261260*u**6+9459450*u**4-18918900*u**2 +4729725)*T**3+(3003*u**10-135135*u**8+1891890*u**6-9459450*u**4+14189175*u**2-2837835)*T**2+(91*u**12-6006*u**10+135135*u**8-1261260*u**6+4729725*u**4-5675670*u**2+945945)*T+u**14-91*u**12+3003*u**10-45045*u**8+315315*u**6-945945*u**4+945945*u**2-135135
  elseif (q.EQ.15) then
     aa= 2027025*u*T**7+(4729725*u**3-14189175*u)*T**6+(2837835*u**5-28378350*u**3+42567525*u)*T**5+(675675*u**7-14189175*u**5+70945875*u**3-70945875*u)*T**4+(75075*u**9-2702700*u**7+28378350*u**5-94594500*u**3+70945875*u)*T**3+(4095*u**11-225225*u**9+4054050*u**7-28378350*u**5+70945875*u**3-42567525*u)*T**2+(105*u**13-8190*u**11+225225*u**9-2702700*u**7+14189175*u**5-28378350*u**3+14189175*u)*T+u**15-105*u**13+4095*u**11-75075*u**9+675675*u**7-2837835*u**5+4729725*u**3-2027025*u 
  elseif (q.EQ.16) then
     aa= 2027025*T**8+(16216200*u**2-16216200)*T**7+(18918900*u**4-113513400*u**2+56756700)*T**6+(7567560*u**6-113513400*u**4+340540200*u**2-113513400)*T**5+(1351350*u**8-37837800*u**6+283783500*u**4-567567000*u**2+141891750)*T**4+(120120*u**10-5405400*u**8+75675600*u**6-378378000*u**4+567567000*u**2-113513400)*T**3+(5460*u**12-360360*u**10+8108100*u**8-75675600*u**6+283783500*u**4-340540200*u**2+56756700)*T**2+(120*u**14-10920*u**12+360360*u**10-5405400*u**8+37837800*u**6-113513400*u**4+113513400*u**2 -16216200)*T+u**16-120*u**14+5460*u**12-120120*u**10+1351350*u**8-7567560*u**6+18918900*u**4-16216200*u**2+2027025
  elseif (q.EQ.17) then
     aa= 34459425*u*T**8+(91891800*u**3-275675400*u)*T**7+(64324260*u**5-643242600*u**3+964863900*u)*T**6+(18378360*u**7-385945560*u**5+1929727800*u**3-1929727800*u)*T**5+(2552550*u**9-91891800*u**7+964863900*u**5-3216213000*u**3+2412159750*u)*T**4+(185640*u**11-10210200*u**9+183783600*u**7-1286485200*u**5+3216213000*u**3-1929727800*u)*T**3+(7140*u**13-556920*u**11+15315300*u**9-183783600*u**7+964863900*u**5-1929727800*u**3+964863900*u)*T**2+(136*u**15-14280*u**13+556920*u**11-10210200*u**9+91891800*u**7-385945560*u**5+643242600*u**3-275675400*u)*T+u**17-136*u**15 +7140*u**13-185640*u**11+2552550*u**9-18378360*u**7+64324260*u**5-91891800*u**3+34459425*u
  elseif (q.EQ.18) then
     aa= 34459425*T**9+(310134825*u**2-310134825)*T**8+(413513100*u**4-2481078600*u**2+1240539300)*T**7+(192972780*u**6-2894591700*u**4+8683775100*u**2-2894591700)*T**6+(41351310*u**8-1157836680*u**6+8683775100*u**4-17367550200*u**2+4341887550)*T**5+(4594590*u**10-206756550*u**8+2894591700*u**6-14472958500*u**4+21709437750*u**2-4341887550)*T**4+(278460*u**12-18378360*u**10+413513100*u**8-3859455600*u**6+14472958500*u**4-17367550200*u**2+2894591700)*T**3+(9180*u**14-835380*u**12+27567540*u**10-413513100*u**8+2894591700*u**6-8683775100*u**4+8683775100*u**2-1240539300)*T**2+(153*u**16-18360*u**14+835380*u**12-18378360*u**10+206756550*u**8-1157836680*u**6+2894591700*u**4-2481078600*u**2+310134825)*T+u**18-153*u**16+9180*u**14-278460*u**12+4594590*u**10- 41351310*u**8+192972780*u**6-413513100*u**4+310134825*u**2-34459425
  elseif (q.EQ.19) then
     aa= 654729075*u*T**9+(1964187225*u**3-5892561675*u)*T**8+(1571349780*u**5-15713497800*u**3+23570246700*u)*T**7+(523783260*u**7-10999448460*u**5+54997242300*u**3-54997242300*u)*T**6+(87297210*u**9-3142699560*u**7+32998345380*u**5-109994484600*u**3+82495863450*u)*T**5+(7936110*u**11-436486050*u**9+7856748900*u**7-  54997242300*u**5+137493105750*u**3-82495863450*u)*T**4+(406980*u**13-31744440*u**11+872972100*u**9-10475665200*u**7+54997242300*u**5-109994484600*u**3+54997242300*u)*T**3+(11628*u**15-1220940*u**13+47616660*u**11-872972100*u**9+7856748900*u**7-32998345380*u**5+54997242300*u**3-23570246700*u)*T**2+(171*u**17-23256*u**15+1220940*u**13-31744440*u**11+436486050*u**9-3142699560*u**7+10999448460*u**5-15713497800*u**3+5892561675*u)*T+u**19-171*u**17+11628*u**15-406980*u**13+7936110*u**11-87297210*u**9+523783260*u**7-1571349780*u**5+1964187225*u**3-654729075*u
  elseif (q.EQ.20) then
     aa= 654729075*T**10+(6547290750*u**2-6547290750)*T**9+(9820936125*u**4-58925616750*u**2+29462808375)*T**8+(5237832600*u**6-78567489000*u**4+235702467000*u**2-78567489000)*T**7+(1309458150*u**8-36664828200*u**6+274986211500*u**4-549972423000*u**2+137493105750)*T**6+(174594420*u**10-7856748900*u**8+109994484600*u**6-549972423000*u**4+824958634500*u**2-164991726900)*T**5+(13226850*u**12-872972100*u**10+19641872250*u**8-183324141000*u**6+687465528750*u**4-824958634500*u**2+137493105750)*T**4+(581400*u**14-52907400*u**12+1745944200*u**10-26189163000*u**8+183324141000*u**6-549972423000*u**4+549972423000*u**2-78567489000)*T**3+ (14535*u**16-1744200*u**14+79361100*u**12-1745944200*u**10+19641872250*u**8-109994484600*u**6+274986211500*u**4-235702467000*u**2+29462808375)*T**2+(190*u**18-29070*u**16+1744200*u**14-52907400*u**12+872972100*u**10-7856748900*u**8+36664828200*u**6-78567489000*u**4+58925616750*u**2-6547290750)*T+u**20-190*u**18+14535*u**16-581400*u**14+13226850*u**12-174594420*u**10+1309458150*u**8-5237832600*u**6+9820936125*u**4-6547290750*u**2+654729075     
  end if
end function aa
